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Abstract 

The Wave Function Matching (WFM) technique has recently been developed for the calculation 
of electronic transport in quantum two-probe systems. In terms of efficiency it is comparable 
with the widely used Green's function approach. The WFM formalism presented so far requires 
the evaluation of all the propagating and evanescent bulk modes of the left and right electrodes 
in order to obtain the correct coupling between device and electrode regions. In this paper we 
will describe a modified WFM approach that allows for the exclusion of the vast majority of the 
evanescent modes in all parts of the calculation. This approach makes it feasible to apply iterative 
techniques to efficiently determine the few required bulk modes, which allows for a significant 
reduction of the computational expense of the WFM method. We illustrate the efficiency of 
the method on a carbon nanotube field-effect-transistor (FET) device displaying band-to-band 
tunneling and modeled within the semi-empirical Extended Hiickel theory (EHT) framework. 

PACS numbers: 73.40.-c, 73.63.-b, 72.10.-d, 85.35.Kt, 85.65. +h 
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FIG. 1: (Color online) Schematic illustration of a nano-scale two-probe system in which a device 
is sandwiched between two semi-infinite bulk electrodes. 

I. INTRODUCTION 



Quantum transport simulations have become an important theoretical tool for investigat- 
ing the electrical properties of nano-scale systems.-'^i^i'^'^ The basis for the approach is the 
Landauer-Biittiker picture of coherent transport, where the electrical properties of a nano- 
scale constriction is described by the transmission coefficients of a number of one-electron 
modes propagating coherently through the constriction. The approach has been used suc- 
cessfully to describe the electrical properties of a wide range of nano-scale systems, including 
atomic wires, molecules and interfaces.-'^'^'^>^»^»^»^>^»^ In order to apply the method to 
semiconductor device simulation, it is necessary to handle systems comprising many thou- 
sand atoms, and this will require new efficient algorithms for calculating the transmission 
coefficient. 

Our main purpose in this paper is to give details of a method we have developed, based on 
the WFM technique,— li^ii^ which is suitable for studying electronic transport in large-scale 
atomic two-probe systems, such as large carbon nanotubes or nano-wire configurations. 

We adopt the many-channel formulation of Landauer and Biittiker to describe electron 
transport in nano-scale two-probe systems composed of a left and a right electrode attached 
to a central device, see Fig. [TJ In this formulation, the conduction Q of incident electrons 
through the device is intuitively given in terms of transmission and reflection matrices, t 
and r, that satisfy the unitarity condition tH + r'''r = 1 in the case of elastic scattering. The 
matrix element tij is the probability amplitude of an incident electron in a mode i in the 
left electrode being scattered into a mode j in the right electrode, and correspondingly 
is the probability of it being reflected back into mode k in the left electrode. This simple 
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interpretation yields the Landauer-Biittiker formula^ 

g = ^Tr[ttt], (1) 

which holds in the limit of infinitesimal voltage bias and zero temperature. 

To our knowledge, the WFM schemes presented so far in the literature require the eval- 
uation of all the Bloch and evanescent bulk modes of the left and right electrodes in order 
to obtain the correct coupling between device and electrode regions. The reason for this is 
that the complete set of bulk modes is needed to be able to represent the proper refiected 
and transmitted wave functions. In this paper we will describe a modified WFM approach 
that allows for the exclusion of the vast majority of the evanescent modes in all parts of the 
calculation. The primary modification can be pictured as a simple extension of the central 
region with a few principal electrode layers. In this manner, it becomes advantageous to 
apply iterative techniques for obtaining the relatively few Bloch modes and slowly decaying 
evanescent modes that are required. We have recently developed such an iterative method 
in Ref. Il9|, which allows for an order of magnitude reduction of the computational expense 
of the WFM method in practice. 

In this work, the proper analysis of the modified WFM approach is presented. The 
accuracy of the method is investigated and appropriate error estimates are developed. As 
an illustration of the applicability of our WFM scheme we consider a 1440 atom CNTFET 
device of 14 nm in length. We calculate the zero-bias transmission curves of the device under 
various gate voltages and reproduce previously established characteristics of band-to-band 
tunneling.-^fl We compare directly the results of the modified WFM method to those of the 
standard WFM method for quantitative verification of the calculations. 

The rest of the paper is organized as follows. The WFM formalism used to obtain t and 
r is introduced in Sect. [Tll In Sect. Hn] we present our method to effectively exclude the 
rapidly decaying evanescent modes from the two-probe transport calculations. Numerical 
results are presented in Sect. [IVl and the paper ends with a short summary and outlook. 

II. FORMALISM 

In this section we give a minimal review of the formalism and notation that is used 
in the current work in order to determine the transmission and refiection matrices t and 
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r. This WFM technique has several attractive features compared to the widely used and 
mathematically equivalent Green's function approach.-!^ Most importantly, the transparent 
Landauer picture of electrons scattering via the central region between Bloch modes of the 
electrodes is retained throughout the calculation. Moreover, WFM allows one to consider the 
significance of each available mode individually in order to achieve more efficient numerical 
procedures to obtain t and r. 



A. Wave function matching 

The WFM method is based upon direct matching of the bulk modes in the left and 
right electrode to the scattering wave function of the central region. For the most part 
this involves two major tasks; obtaining the bulk electrode modes and solving a system of 
linear equations. The bulk electrode modes can be characterized as either propagating or 
evanescent (exponentially decaying) modes but only the propagating modes contribute to Q 
in Eq. ([T]). We may write Q = {2e'^/h)T, where 

T = J2\tkk'\' (2) 

kk' 

is the total transmission and the sum is limited to propagating modes k and k' in the left and 
right electrode, respectively. Notice, however, that the evanescent modes are still needed in 
order to obtain the correct matrix elements tkk'- We will discuss this matter in Sect. IIII CI 
We assume a tight-binding setup for the two-probe systems in which the infinite structure 
is divided into principal layers numbered i = — oo, . . . , oo and composed of a finite central 
(C) region containing the device and two semi-infinite left (L) and right (R) electrode 
regions, see Fig. [2j The wave function is il^i{x) = CijXiji^ ~ Xjj) in layer i, where Xi,j 
denotes localized non-orthogonal atomic orbitals and Xj ,,• are the positions of the rrii orbitals 
in layer i. We represent ij^i{x) by a column vector of the expansion coefficients, given by 
V'i = • ■ ■ ,Ci^miV, and write the wave function V extending over the entire system as 
■^A = [V'-oo) • • • 5 We also assume that the border layers 1 and n of the central region 



are always identical to a layer of 
We refer the reader to Refs 



;he connecting electrodes. 

M to our setup. 

Here and in the rest of this paper, we will use the following notation for the key elements: 
The matrices = . . . , ^l,™^] contain in their columns the full set of rriL left-going 
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FIG. 2: (Color online) Schematic representation of WFM applied to layered two-probe systems, 
where the central device region, consisting of layers i = l,...,n, is attached to left and right 
semi-infinite electrodes. The incoming propagating mode from the left electrode is scattered in 
the central region and ends up as reflected and transmitted superpositions of propagating and 
evanescent modes. 

(— ) and rriL right-going (-|-) bulk modes (pf^j^ of the left electrode, and the diagonal matrices 
= diag[A^^, X^^^ ■ ■ ■ ^ -^l mi] ^o\d the corresponding Bloch factors. If trivial modes with 
\^Lk\ = or = C50 occur they are simply rejected. We assume that all the evanescent 

bulk modes are (state-)normalized (f'^^fp'tk ~ '^hile all the Bloch bulk modes are flux- 
normalized^ 4>i,\(f}i,k ~ ^L/v'^k^ where f are the group velocities'"^ and di is the layer 
thickness. Similarly for the right electrode the matrices and are formed. 

We also introduce the Bloch matrices^^ = $^A^($^)~^ and = $^A^($^)~^ 
which propagate the layer wave functions in the bulk electrode 

= (B±)^-#, (3) 

where subscript L is implied for the left electrode < 1), and R for the right electrode 
> n). Notice that the first central region layer is defined for layer 1 and not layer 0, as 
is the case in Ref. 3 

As explicitly shown in Refs. 1^,17 isl, by fixing the layer wave functions coming into 



the C region (e.g., in our case i/^i = Xj^k^L^ '^n = 0) and matching the layer wave 
functions across the C region boundaries, the system of linear equations for the central 
region wavefunction il^c can be written as 

{ESc - He - Si - 1:r)^c = b, (4) 
where E is the energy, the overlap and the Hamiltonian matrix of the central region. 



In the following we discuss the terms, I]^, "Er, and b, which arise from matching the 
boundary conditions with the electrode modes. 

The self-energy matrices, S ^ and S/j, arise from matching with the outgoing left and right 
electrode modes. They only have non-zero terms in the upper left and lower right corner 
block, respectively, and these elements can be calculated in terms of the Bloch matrices:— 

[Si]i,i = hS,i(Hi + fit i(B^)-i)-iHo,i, (5) 

and 

(6) 

where we have introduced the overline notation fij = ESi — Hj and fijj- = ESij — Hij- For 
the current setup, these matrices are identical to the self-energy matrices introduced in the 
Green's function formalism^ (to within an infinitesimal imaginary shift of E), and may be 
evaluated by well-known recursive techniques^^i^i or constructed directly from the electrode 
modes using Eq. (Q. 

The source term b arises from the incoming mode. Assuming an incoming mode from 
the left, we have b = [hj, O"^, . . . , O"""]""" specified by the expression 

bi = -{Ul, + [Si]i,iB+)V>o, (7) 

where i/^o is the incoming wave function. 

For notational simplicity in the following sections, we leave out the implied subscripts 
L or R, indicating the left or right electrode, whenever the formalism is the same for both 
(e.g, for symbols m, Xk, 4>k, A^, B^, S, etc.). 

B. Transmission and reflection coefficients 

As a final step we want to determine the t and r matrices from the boundary wave 
functions 'j/'i and i/^n that have been obtained by solving Eq. (jlj). 

When the incoming wave Vo is specified to be the kth right-going mode (p'j^^. of the left 
electrode, then il^n will be the superposition of outgoing right transmitted waves. The fcth 
column of the transmission matrix is defined as the corresponding expansion coefficients 
in right electrode modes and can be evaluated by solving 

$^tfc = (8) 
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TABLE I: CPU times in seconds when using WFM for calculating t and r at 20 different energies 
inside E G [—2 eV; 2 eV] for various two-probe systems. The numbers of atoms in the central 
region (electrode unit cell) are indicated. The four right-most columns show the CPU times spent 
for computing the electrode bulk modes with dgeev and in this work vs. solving the central region 
linear systems in Eq. ^ and the system with two extra principal layers on each side. 
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Al-Cx7-Al 


74(18) 


0.4 


0.6 


3.6 


1.6 


Au-DTB-Au 


102(27) 


8.1 


13.5 


91.0 


28.2 


Au-CNT(8,0)xl-Au 


140(27) 


11.4 


16.6 


77.6 


17.1 


Au-CNT(8,0)x5-Au 


268(27) 


45.3 


50.3 


83.6 


17.8 


CNT(8,0)-CNT(8,0) 


192(64) 


7.0 


11.9 


129.0 


19.4 


CNT(4,4)-CNT(8,0) 


256(64|64) 


7.2 


12.4 


121.5 


21.0 


CNT(5,0)-CNT(10,0) 


300(40|80) 


24.7 


31.5 


113.3 


22.6 


CNT(18,0)-CNT(18,0) 


576(144) 


172.2 


225.5 


1362.2 


253.3 


CNTFET (see Fig. E]) 


1440(160) 


259.8 


286.9 


4633.0 


372.3 



where $J is the rriR x m/j column matrix holding the right-going bulk modes of the right 
electrode (and here assumed to be non- singular). Similarly the kth column of the reflection 
matrix is given by 

$^rfc = V'l - A+,0+„ (9) 

where holds the left-going bulk modes of the left electrode. The flux normalization 
ensures that tH + r^r = 1. 

III. EXCLUDING EVANESCENT MODES 

The most time consuming task of the WFM method is often to determine the electrode 
modes, which requires solving a quadratic eigenvalue problem.-i^ As examples, see the pro- 
filing results listed in Table [H where we have used the method to compute t and r for a 
selection of two-probe systems.— The CPU timings show that to determine the electrode 
modes by employing the state-of-the-art lapack eigensolver dgeev is, in general, much 
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more expensive than to solve the system of hnear equations in Eq. (jlj). We expect this 
trend to hold for larger systems as well. Therefore, in the attempt to model significantly 
larger devices (thousands of atoms), it is of essential interest to reduce the numerical cost 
of the electrode modes calculation. We argue that a computationally reasonable approach 
is to limit the number of electrode modes taken into account, e.g., by excluding the least 
important evanescent modes. In this section, a proper technique to do this in a rigorous and 
systematic fashion is presented. 

A. Decay of evanescent modes 

The procedure to determine the Bloch factors and non-trivial modes 4>k of an ideal 
electrode and subsequently characterize these as right-going (+) or left-going (— ) is well 
described in the literature.— ii^>i^>2^ We note that only the obtained propagating modes with 
\Xk\ = 1 are able to carry charge deeply into the electrodes and thus enter the Landauer 
expression in Eq. ([2]). The evanescent modes with |Afc| 7^ 1, on the other hand, decay 
exponentially but can still contribute to the current in a two-probe system, as the "tails" 
may reach across the central region boundaries. 

Consider a typical example of an electrode modes evaluation: We look at a gold electrode 
with 27 atoms in the unit cell represented by 9 (sp^d^) orbitals for each Au-atom. Such a 
system results in 243 right-going and 243 left-going modes. Fig. [3K shows the positions in the 
complex plane of the Bloch factors corresponding to the right-going modes (i.e., |Afc| < 1) 
for energy E = —1.5 eV. We see that there are exactly three propagating modes, which 
have Bloch factors located on the unit circle. The remaining modes are evanescent, of which 
many have Bloch factors with small magnitude very close to the origin. 

Fig. [2)3 illustrates how the 243 left-going modes would propagate through 10 successive 
gold electrode unit cells. The figure shows that the amplitudes of the three propagating 
modes are unchanged, while the evanescent modes are decaying exponentially. In particular, 
we note that the evanescent modes with Bloch factors of small magnitude are very rapidly 
decaying and vanishes in comparison to the propagating modes after only a few layers. 
In the following, we will exploit this observation and attempt to exclude such evanescent 
modes from the WFM calculation altogether. Formally this can be accomplished if only the 
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FIG. 3: (Color online) (a) Positions of the Bloch factors (|Afc| < 1) obtained for a bulk Au(lll) 
electrode with 27 atoms per unit cell at E = —1.5 eV. (b) Amplitudes of the corresponding 
normalized electrode modes cj)k moving through 10 layers of the ideal bulk electrode. A total of 
243 modes are shown of which 3 are propagating (colored/dashed) and the rest are evanescent 
(circles/black). 

electrode modes 0^ with Bloch factors satisfying 

Amin < l^fcl < -^min' (lO) 

are computed and subsequently taken into account, for a reasonable choice of < Amin < 1- 
Eq. (fTOj) is adopted as the key relation to select a particular subset of the available electrode 



modes (as recently suggested in Ref. llTl ) 



B. Extra electrode layers 

We will denote the mode, Bloch and self-energy matrices from which the rapidly decaying 
evanescent modes are excluded with a tilde, i.e., as B='= and S. The mode matrices 
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FIG. 4: (Color online) Two-probe system in which the C region boundaries are expanded by I 
extra electrode layers. 

holding the excluded modes are denoted by a math-ring accent so that 

$±=[$±,$±]^ (11) 

is the assumed splitting of the full set. All expressions to evaluate the Bloch and self-energy 
matrices are unchanged as given in Sect. [TTl (now ($^)~^ merely represents the pseudo- 
inverses of $^). However, since the column spaces of are not complete, there is no 
longer any guaranty that WFM can be performed so that the resulting self-energy matrices 
and, in turn, the solution ipc = [V'?) • • • ; of ^^e linear system in Eq. (jlj), are correct. 
In addition, it is clear that errors can occur in the calculation of t and r from Eqs. ([8]) and 
([9]) because the boundary wave functions i/'i and i/'n might not be fully represented in the 
reduced sets $J and 

In order to diminish the errors introduced by excluding evanescent modes, we propose 
to insert additional electrode layers in the central region, see Fig. |H As illustrated in the 
previous section, this would quickly reduce the imprint of the rapidly decaying evanescent 
modes in the boundary layer wave functions Vi ^^^d i/'n) which means that the critical 
components outside the column spaces becomes negligible at an exponential rate in terms 
of the number of additional layers. We emphasize that the inserted layers may be "fictitious" 
in the sense that they can be accommodated by simple block-Gaussian-eliminations prior to 
the solving of Eq. (jlj) for the original system. 

The above statements are confirmed by the following analysis. We expand the electrode 
wave functions in the corresponding complete set of bulk modes 



(12) 
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where = [a^^, a^^]""" are vectors that contain the expansion coefficients. In the particular 
case, where I extra electrode layers are inserted and the border layers of the C region 
are identical to the connecting electrode layers, the electrode wavefunctions entering the 
matching boundary equations will be 



(0- 



(A 
(A, 



7)-'a7 



(13) 



and 



(14) 



using the definition B 



± 



(A^3'a+ 
^ (A^^)'a+ 

$^A^($^)^^. This shows that the critical components out 



side the column spaces of and are given by coefficients (Az)~'aj~ and (A^)'a^, 
respectively. If this set only consists of the most rapidly decaying of the evanescent modes 
according to Eq. (fTOl) . that is, |Aa,.| > A^j^ for the diagonal elements of and |Aa,.| < Amin 
for the diagonal elements of A^, where Amin is less than 1, these coefficients always decrease 
as a function of /. 

We conclude that WFM with the reduced set of modes approaches the exact case if 
additional electrode layers are inserted and the solution i/'c obtained from Eq. (jl]) approaches 
the correct solution il^c accordingly. 



C. Accuracy 

As pointed out above, the exclusion of some of the evanescent modes from the mode 
matrices will introduce errors because the column spaces in are incomplete. In this 
section we will estimate how this will infiuence the accuracy of the calculated transmission 
and refiection coefficients in terms of the parameter Amin and the number / of extra electrode 
layers. 

Consider first the accuracy of the transmission matrix t in the case of the extended two- 
probe system in Fig. |H For a specific incoming mode k, we compare the correct result 
obtained with the complete set of modes (cf. Eq. (IHl)), 



tfc 



[^^R,^^^]-'i^^^\ (15) 



11 



with the result obtained with the reduced mode matrix (denoted by a prime), 



$+6]-iVi')+, (16) 



0' 

where 0' represents the zero vector of size and the zero matrix of size mn x m^. 

The important coefficients in and t'^ for transmission calculations are the ones repre- 
senting the Bloch modes which enters the Landauer-Biittiker formula in Eq. ([2]). Since these 
are never excluded they will always be located within the first elements, i.e., in and 
t'f.. It then suffices to compare these parts of the transmission matrix which we can do as 
follows. 

From the properties of the pseudo-inverse we are able to write the relation 

{Hr'[^tMn] = [i{Hr'^ti (17) 

where I is the identity matrix of order equal to the number of included modes ttlr- Using 
the expression in Eq. f|T4|) it then follows that 

tfc = (A+)'a+, (18) 

and 

t', = t, + {^V)-'H{k-^n)'K. (19) 
where the expression clearly corresponds to the correct coefficients plus an error term. 

We have already established in the previous section that the [Al^^^k^ factor in the error 
term will decrease as a function of /. We now show that the other term, ($^)~^$^ is 
independent of Z, and consequently, that the error term in Eq. f[T^ must decrease as a 
function of I. To this end we look at the 2-norm of ($^)~^$^, which satisfies 

Il(^^^)-'^^il|2<4ll(^^)-^ll2, (20) 

□ 1_ 

since ||$^||2 < when all evanescent modes are assumed to be normalized. The norm 
||($^)^^||2 can be readily evaluated and depends on the set of modes included via the 
parameter Amin but not on /. Thus, we conclude that the only term of Eq. f|T9l) which 
depend on / is (A^)'a^, and the error is therefore decreasing as function of I. 

Writing Eq. ( fT9|) as tj, = + e^, where holds the errors on the coefficients of the fcth 
column, we further obtain that the total transmission T' can be expressed as 

T' = T + ^(t^fc,efcfc' + ^Ik'^kk' + lefcfe'H (21) 



kk' 
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where T is the exact result and the summation is over the Bloch modes k and k' in the left 
and right electrode, respectively. 

For a first order estimate of the error term in Eq. (1211) we consider the worst case ap- 
proximation, where all diagonal elements of are equal to the maximum range Amin of 
Eq. fllOl) . This makes all elements e^k' proportional to Aj^j^. and we arrive at the simple 
relation 

\T'-T\^Xl,, + 0{{XU'), (22) 

which shows that the error decreases exponentially in terms of the number of extra layers I. 

For a higher order estimate of the error, we directly monitor the error arising on the 
boundary conditions, in terms of the coefficient vectors = {^r)~^{'^?^ — ^tk^'t k) ^^"^ 
bR,k = {^R)~^'4^n''~ , where and i/'n^" are given by solving Eq. (jl]). When the boundary 
conditions are exactly satisfied, we have 1^^^! = and |b_R,fc| = 0. In the case where the 
boundary conditions are not exactly satisfied, h^^k represents the error on the left-going 
components within the right boundary layer in the same way that represents the error 
on the right-going (transmitted) components. We would therefore expect the same order of 
magnitude of |b^fc| and \ek\ in an actual calculation for a given mode k. This suggests the 
following error estimate from Eq. ( !2TI) . 

\r -T\< ^(2|tfc||efc| + 141') ~ J2^2\ik\\bR^k\ + \bR,k\^), (23) 

k k 

where all the vector norms (e.g., jt^p = ^j^, It^^'p) are assumed to be taken over the 
elements corresponding to Bloch bulk modes k' only. 

Finally, we note without explicit derivation, that similar arguments for the reflection 
matrix with columns r'f^ = ($^)^^(i/'f^ — Xj^/^fptk) total reflection coefficient R', 

results in the same accuracy expressions for \R' — R\ if we substitute and bR^^ ^L,fe 

in Eqs. (j22|) and (|23|). 

D. Example 

To end this section, we exemplify the previous discussion quantitatively by looking at the 
Au(lll) electrode described earlier, and assuming a 128 atom (4 unit cells) device of zigzag- 
(8,0) carbon nano tube (CNT) sandwiched between the gold electrodes, see the configuration 
in Fig. [H For energy E = —1.5 eV, we have calculated the deviation between the total 
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FIG. 5: (Color online) Error (absolute) in the calculated total transmission (solid red lines) and 
reflection (solid blue lines) coefficients T' and R' as a function of /. The panels show the cases of 
Amin set to 0.5, 0.3 and 0.1, which corresponds to 3, 14 and 31 Au bulk modes (out of 243, see 
Fig. [3]) taken into account, respectively. The dashed line indicates the first order error estimate 
'^min- yellow and green lines show error estimates obtained from Eq. (j23j) . 

transmission obtained when all bulk modes are taken into account (T) and when some 
evanescent modes are excluded (T') as specified with different settings of Amin- Deviations 
are also determined for the corresponding total reflection coefficients {R and R'). Fig. [5] 
shows the results as a function of Z, together with the estimate A^j^ of Eq. fl22|) and the 
estimate of Eq. (1231) both for the transmission and reflection coefficients, where the higher 
order terms have been neglected. 

We observe that the absolute error in the obtained transmission coefficients (red curves) 
and reflection coefficients (blue curves) are generally decreasing as a function of /, following 
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the same convergence rate as Aj^j^ (dashed hne). Looking closer at results for neighbor / 
values, we see that the errors initially exhibit wave-like oscillations. This is directly related 
to the wave form of the evanescent modes that have been excluded (see the propagation of 
the slowest decaying black curves in Fig. El^b)). In other words, although the norm of the 
errors |efc| are decreasing as a function of /, the specific error ikk' on a given (large) coefficient 
of t'l^i^, or f^^, may increase, which means that the overall error term in Eq. (!2T!) can go up. 
Fortunately this is only a local phenomenon with the global trend being rapidly decreasing 
errors. 

Consider also the quality of the simple accuracy estimate of Aj^jn and the estimates ex- 
pressed by Eq. f l23l) for the transmission coefficients (green curves) and reflection coefficients 
(yellow curves), respectively. For relatively large Amin all estimates are very good. However, 
for smaller values of Amin, only the latter two retain a high quality while the Aj^j^ estimate 
tends to be overly pessimistic. It is important to remember that these estimates are by no 
means strict conditions but in practice give very reasonable estimates of the accuracy. 

We note in passing, that the results in the top panel of Fig. [5] corresponds to using only 
the propagating Bloch modes in the transmission calculation. Still we are able to compute 
T and R to an absolute accuracy of three digits by inserting 2x5 extra electrode layers in 
the two-probe system. This is quite remarkable and shows promise for large-scale systems, 
e.g., with nano-wire electrodes, for which the total number of evanescent modes available 
becomes exceedingly great. 

IV. APPLICATION 

In this section we will apply the developed method to a nano-device consisting of a CNT 
stretched between to two metal electrodes and controlled by three gates. The setup is 
inspired by Appenzeller et al.,^ and we expect this particular arrangement to be able to 
display so-called band-to-band (BTB) tunneling, where one observes gate induced tunneling 
from the valence band into the conduction band of a semi-conducting CNT and vice versa. 

We show the conflguration of the two-probe system in Fig. [61 The device conflguration 
contains 10 principal layers of a CNT(8,4), having 112 atoms in each layer. The diameter 
of the tube and the thickness of the principal layer are 8.3 A and 11.3 A, respectively. The 
electrodes consist of CNT(8,4) resting on a thin surfaces of Li, where the lattice constant 
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FIG. 6: (Color online) Schematic illustration of a carbon nanotube (8,4) band-to-band tunneling 
device. The carbon nanotube is positioned on Li surfaces next to an arrangement of three gates. 

of the Li layers is stretched to fit the layer thickness of the CNT. The central region of the 
two-probe system comprises a total of 1440 atoms. An arrangement of rectangular gates 
are positioned below the carbon nanotube as indicated on the figure. In the plane of the 
illustration (length x height) the dimensions are as follows: Dielectric 108 A x 5 A; Gate-A 
108 A X 5 A; Gate-B 20 A x 5 A. We set e = 4 for the dielectric constant of the dielectric 
in order to simulate Si02 or AI2O3 oxides. All the regions are centered with respect to 
the electrodes so that the complete setup has mirror symmetry in the length direction. In 
the direction perpendicular to the illustration the configuration is assumed repeated every 
19.5 A as a super-cell. 

We have obtained the density matrix of the BTB device by combining the NEGF for- 
malism with a semi-empirical Extended Hiickel model (EHT) using the parameterization 
of Hoffmann. From the density matrix we calculate Mulliken populations on each atom, 
and represent the total density of the system as a superposition of Gaussian distributions 
on each atom properly weighted by the Mulliken population. The width of the Gaussian is 
chosen to be consistent with CNDO parameters.— The electrostatic interaction between the 
charge distribution and the dielectrics and gates is subsequently calculated. The Hartree- 
like term is then included in the Hamiltonian and the combined set of equations are solved 
self-consistently. The resulting self-consistent EHT model is closely related to the work of 



Ref. l27|, and a detailed description of the model will be presented elsewhere.— 

In order to adjust the charge transfer between the CNT and the Li electrodes we add 
the term 5eS to the Li parameters. With an appropriate adjusted value of (5e, the carbon 
nanotube becomes n-type doped. We adjust the value such that the average charge transfer 
from Li to the nanotube at self-consistency is 0.002 e per carbon atom in the electrode. The 
Fermi energy is then located at —4.29 eV, which is 0.07 eV below the conduction band of 
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FIG. 7: (Color online) Left panel: Representation of the electrostatic induced shift of the valence 
and conduction band edges along the length of the device for gate potentials Vcatc-B = —2.0 eV, 
1.0 eV, 2.0 eV and 4.0 eV. Right panel: The corresponding transmission spectrum. The dotted 
line shows the position of the Fermi level, and the solid line shows the transmission coefficient for 
an ideal CNT(8,4). 

the CNT(8,4). 

In the following we fix Vcatc-A = —2.0 eV and vary the Gate-B potentials in the range 
[—2 eV, 4 eV]. Note that we report the gate potentials as an external potential on the 
electrons, and to translate the values into a gate potential of unit Volts the values must be 
divided with — e. 

In the left part of Fig. [7] we present the total self-consistent potential induced by the 
three gates on the carbon atoms in the CNT over the full extension of the device. For each 
configuration of the gate potentials the electrostatic potential is shown twice, i.e., by two 
curves with the same color displaced relative to each other with the energy of the valence 
band and conduction band edge, respectively. In this way the curves not only represent 
the electrostatic potential of the device, but also the position of the valence and conduction 
band edges. 

Along with this, in the right part of Fig. [71 we show the corresponding transmission 
spectrum T{E), for four gate potentials Vcate-B = —2.0 eV, 1.0 eV, 2.0 eV, and 4.0 eV. When 
Vcate-B = —2.0 eV the nanotube is largely unpertubed by the gate and the transmission 
coefficient is close to an ideal (8,4) CNT. We note that this is in agreement with ab initio 
calculations by Nardelli et. al.,^^ which found that a two terminal (5,5) CNT device in 
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FIG. 8: (Color online) Conduction in units of the conductance quantum Go as a function of the 
Gate-B potential. In the calculations we use a dielectric constant of 4, Vcate-A = —2.0 eV, and 
vary Vcate-B from —2.0 eV to 4.0 eV as indicated. 

a similar contact geometry showed a nearly ideal conductance spectrum. In addition, the 
calculated band gap of the (8,4) nanotube is 0.81 eV, which is in good agreement with the 
value of 0.96 eV obtained from ab initio density-functional calculations in the generalized 
gradient approximation.— 

From Fig. [7] we see how the bands are shifted upwards by an increasing amount as the 
Gate-B potential is turned up. To begin with, e.g., for Vcate-B = 1 eV, this results in 
lower conduction since the conduction band bends away from the Fermi level and the Fermi 
energy electrons need to tunnel through the central region. When the gate voltage is at 
^Gate-B = 2 cV, the valcucc band almost reaches the conduction band in which case BTB 
tunneling becomes possible. By increasing the gate voltage further, more bands become 
available for BTB tunneling and the effect is visible as a steady increase in the calculated 
transmission T{E) just above the Fermi level. 

The results for the Fermi level transmission T{Ef) corresponding to the T = K unit 
conduction Go, are displayed with the black curve in Fig. [HI It shows an initial conductance 
for Vcatc-B = —2.0 V of the order of one, a subsequent drop by four orders of magnitude 
around Vcate-B = 2.0 V, and a final increase of one order of magnitude towards Vcate-B = 
4.0 V. We also display the results for the room temperature T = 300 K conductance (red 
curve), which can be obtained from 



The two conduction curves are similar, showing that the device is operating in the tunneling 




(24) 
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FIG. 9: (Color online) Transmission coefficients T and T calculated with the standard WFM 
method (black solid) and the method of this work (red dashed), respectively, and the difference 
|T — T\ (blue line) as a function of energy E in the Vcate-B = 2 V case. 

regime rather than the thermal emission regime. 

We next briefly comment on the comparison of the simulation with the experiment of 
Appenzeller et al.M^ In both cases the conduction curves have two branches, which we denote 
Field Emission (FE) and Band to Band Tunneling (BTB). Initially, the conduction decreases 
with applied gate potential due to the formation of a barrier in the central region, this is the 
FE regime. For larger biases the conduction increases again due to BTB tunneling, this is 
the BTB regime. The experimental device display thermal emission conduction and shows 
a corresponding subthreshold slope, S, of kBTln{10)/e ~ 60 mV/dec in the FE regime. 
The theoretical device, on the other hand, display tunneling conduction and has S ~ 500 
mV/dec in the FE regime. In the BTB regime, the theoretical device has 5* ~ 2000 mV/dec, 
while the experimental device show S ^ 40 mV/dec. 

The very different behavior is due to the short channel length of the theoretical device. 
The central barrier has a length of ~ 5 nm, and at this length the electron can still tunnel 
through the barrier. We see that the short channel length not only affects the subthreshold 
slope of the FE regime, but also strongly influence the BTB regime. Work are in progress 
for a parallel implementation of the methodology, which will make it feasible to simulate 
larger systems, and thereby investigate the transition from the tunnelling to the thermal 
emission regime. 

All the above results have been calculated with the modified WFM method using pa- 



19 



rameters Amin = 0.1 and 1 = 1. Thus, the results presents a non-trivial application of the 
new method. To verify the transmission results in Fig. [7| we present a comparison with 
the standard WFM method in Fig. O The figure shows that the transmissions curves are 
identical to about three significant digits. The CPU time required for calculating a com- 
plete transmission spectrum for Fig. [7] is (~ 3 hours), while the corresponding calculation 
presented in Fig. [9] with the standard WFM method took (~ 35 hours). Thus, the overall 
time saving achieved with the new method was therefore more than an order of magnitude. 
The results in Table |T] indicate that similar timesavings can be expected for other systems 
with non-trivial electrodes. 



V. SUMMARY 

We have developed an efficient approach for calculating quantum transport in nano- 
scale systems based on the WFM scheme originally proposed by Ando in reference [l6j. In 
the standard implementation of the WFM method for two-probe systems, all bulk modes 
of the electrodes are required in order to represent the transmitted and reflected waves 
in a complete basis. By extending the central region of the two-probe system with extra 
electrode principal layers, we are able to exclude the vast majority of the evanescent bulk 
modes from the calculation altogether. Our final algorithm is therefore highly efficient, and 
most importantly, errors and accuracy can be closely monitored. 

We have applied the developed WFM algorithm to a CNTFET in order to study the 
mechanisms of band-to-band tunneling. The setup was inspired by reference 20|] , and the 
calculation display features also observed in the experiment, however, due to the short chan- 
nel length the theoretical device operates in the tunneling regime, while the experimental 
device operates in the thermal emission regime. 

By measuring the CPU-times for calculating transmission spectra of the CNTFET two- 
probe system and comparing to the cost of the standard WFM method we have observed a 
speed-up by more than a factor of 10. We see similar speed-up for other non-trivial systems. 
We therefore believe that this is an ideal method to be used with ab-initio transport schemes 
for large-scale simulations. 
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Bloch's theorem^-^ xpi = Afci/^j-i for the ideal electrodes defines the phase factors = e"'*'^, 
where is the complex wave number and d is the layer thickness, which are referred to as 
Bloch factors throughout this paper. 

When using the Landauer formula in Eq. ([1]) it is assumed that the electrode Bloch modes 
carry unit current in the conduction direction. This can be conveniently accommodated by flux- 
normalizing the Bloch modes, i.e., <^^^ — > (d-L / jj^ '^^ ™ '^^^ electrode.— 
We should point out that the metallic electrodes in the two-probe systems considered in Table U] 
can be fully described by much smaller unit cells than indicated (often only a few atoms are 
needed) and therefore the time spend on computing the bulk modes can be vastly reduced in 
these specific cases. For a general method, however, which supports CNTs, nano wires, etc. as 
electrodes, the timings are appropriate for showing the overall trend in the computational costs. 
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